home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 1 Issue 2 / PDCD-1 - Issue 02.iso / _utilities / utilities / 001 / meschach / !Meschach / c / zvecop < prev   
Text File  |  1994-03-08  |  11KB  |  511 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26.  
  27. #include    <stdio.h>
  28. #include    "matrix.h"
  29. #include    "zmatrix.h"
  30. static    char    rcsid[] = "$Id: zvecop.c,v 1.2 1994/03/08 05:51:07 des Exp $";
  31.  
  32.  
  33.  
  34. /* _zin_prod -- inner product of two vectors from i0 downwards
  35.     -- flag != 0 means compute sum_i a[i]*.b[i];
  36.     -- flag == 0 means compute sum_i a[i].b[i] */
  37. complex    _zin_prod(a,b,i0,flag)
  38. ZVEC    *a,*b;
  39. u_int    i0, flag;
  40. {
  41.     u_int    limit;
  42.  
  43.     if ( a==ZVNULL || b==ZVNULL )
  44.         error(E_NULL,"_zin_prod");
  45.     limit = min(a->dim,b->dim);
  46.     if ( i0 > limit )
  47.         error(E_BOUNDS,"_zin_prod");
  48.  
  49.     return __zip__(&(a->ve[i0]),&(b->ve[i0]),(int)(limit-i0),flag);
  50. }
  51.  
  52. /* zv_mlt -- scalar-vector multiply -- may be in-situ */
  53. ZVEC    *zv_mlt(scalar,vector,out)
  54. complex    scalar;
  55. ZVEC    *vector,*out;
  56. {
  57.     /* u_int    dim, i; */
  58.     /* complex    *out_ve, *vec_ve; */
  59.  
  60.     if ( vector==ZVNULL )
  61.         error(E_NULL,"zv_mlt");
  62.     if ( out==ZVNULL || out->dim != vector->dim )
  63.         out = zv_resize(out,vector->dim);
  64.     if ( scalar.re == 0.0 && scalar.im == 0.0 )
  65.         return zv_zero(out);
  66.     if ( scalar.re == 1.0 && scalar.im == 0.0 )
  67.         return zv_copy(vector,out);
  68.  
  69.     __zmlt__(vector->ve,scalar,out->ve,(int)(vector->dim));
  70.  
  71.     return (out);
  72. }
  73.  
  74. /* zv_add -- vector addition -- may be in-situ */
  75. ZVEC    *zv_add(vec1,vec2,out)
  76. ZVEC    *vec1,*vec2,*out;
  77. {
  78.     u_int    dim;
  79.  
  80.     if ( vec1==ZVNULL || vec2==ZVNULL )
  81.         error(E_NULL,"zv_add");
  82.     if ( vec1->dim != vec2->dim )
  83.         error(E_SIZES,"zv_add");
  84.     if ( out==ZVNULL || out->dim != vec1->dim )
  85.         out = zv_resize(out,vec1->dim);
  86.     dim = vec1->dim;
  87.     __zadd__(vec1->ve,vec2->ve,out->ve,(int)dim);
  88.  
  89.     return (out);
  90. }
  91.  
  92. /* zv_mltadd -- scalar/vector multiplication and addition
  93.         -- out = v1 + scale.v2        */
  94. ZVEC    *zv_mltadd(v1,v2,scale,out)
  95. ZVEC    *v1,*v2,*out;
  96. complex    scale;
  97. {
  98.     /* register u_int    dim, i; */
  99.     /* complex    *out_ve, *v1_ve, *v2_ve; */
  100.  
  101.     if ( v1==ZVNULL || v2==ZVNULL )
  102.         error(E_NULL,"zv_mltadd");
  103.     if ( v1->dim != v2->dim )
  104.         error(E_SIZES,"zv_mltadd");
  105.     if ( scale.re == 0.0 && scale.im == 0.0 )
  106.         return zv_copy(v1,out);
  107.     if ( scale.re == 1.0 && scale.im == 0.0 )
  108.         return zv_add(v1,v2,out);
  109.  
  110.     if ( v2 != out )
  111.     {
  112.         tracecatch(out = zv_copy(v1,out),"zv_mltadd");
  113.  
  114.         /* dim = v1->dim; */
  115.         __zmltadd__(out->ve,v2->ve,scale,(int)(v1->dim),0);
  116.     }
  117.     else
  118.     {
  119.         tracecatch(out = zv_mlt(scale,v2,out),"zv_mltadd");
  120.         out = zv_add(v1,out,out);
  121.     }
  122.  
  123.     return (out);
  124. }
  125.  
  126. /* zv_sub -- vector subtraction -- may be in-situ */
  127. ZVEC    *zv_sub(vec1,vec2,out)
  128. ZVEC    *vec1,*vec2,*out;
  129. {
  130.     /* u_int    i, dim; */
  131.     /* complex    *out_ve, *vec1_ve, *vec2_ve; */
  132.  
  133.     if ( vec1==ZVNULL || vec2==ZVNULL )
  134.         error(E_NULL,"zv_sub");
  135.     if ( vec1->dim != vec2->dim )
  136.         error(E_SIZES,"zv_sub");
  137.     if ( out==ZVNULL || out->dim != vec1->dim )
  138.         out = zv_resize(out,vec1->dim);
  139.  
  140.     __zsub__(vec1->ve,vec2->ve,out->ve,(int)(vec1->dim));
  141.  
  142.     return (out);
  143. }
  144.  
  145. /* zv_map -- maps function f over components of x: out[i] = f(x[i])
  146.     -- _zv_map sets out[i] = f(x[i],params) */
  147. ZVEC    *zv_map(f,x,out)
  148. #ifdef PROTOYPES_IN_STRUCT
  149. complex    (*f)(complex);
  150. #else
  151. complex (*f)();
  152. #endif
  153. ZVEC    *x, *out;
  154. {
  155.     complex    *x_ve, *out_ve;
  156.     int    i, dim;
  157.  
  158.     if ( ! x || ! f )
  159.         error(E_NULL,"zv_map");
  160.     if ( ! out || out->dim != x->dim )
  161.         out = zv_resize(out,x->dim);
  162.  
  163.     dim = x->dim;    x_ve = x->ve;    out_ve = out->ve;
  164.     for ( i = 0; i < dim; i++ )
  165.         out_ve[i] = (*f)(x_ve[i]);
  166.  
  167.     return out;
  168. }
  169.  
  170. ZVEC    *_zv_map(f,params,x,out)
  171. #ifdef PROTOTYPES_IN_STRUCT
  172. complex    (*f)(void *,complex);
  173. #else
  174. complex    (*f)();
  175. #endif
  176. ZVEC    *x, *out;
  177. void    *params;
  178. {
  179.     complex    *x_ve, *out_ve;
  180.     int    i, dim;
  181.  
  182.     if ( ! x || ! f )
  183.         error(E_NULL,"_zv_map");
  184.     if ( ! out || out->dim != x->dim )
  185.         out = zv_resize(out,x->dim);
  186.  
  187.     dim = x->dim;    x_ve = x->ve;    out_ve = out->ve;
  188.     for ( i = 0; i < dim; i++ )
  189.         out_ve[i] = (*f)(params,x_ve[i]);
  190.  
  191.     return out;
  192. }
  193.  
  194. /* zv_lincomb -- returns sum_i a[i].v[i], a[i] real, v[i] vectors */
  195. ZVEC    *zv_lincomb(n,v,a,out)
  196. int    n;    /* number of a's and v's */
  197. complex    a[];
  198. ZVEC    *v[], *out;
  199. {
  200.     int    i;
  201.  
  202.     if ( ! a || ! v )
  203.         error(E_NULL,"zv_lincomb");
  204.     if ( n <= 0 )
  205.         return ZVNULL;
  206.  
  207.     for ( i = 1; i < n; i++ )
  208.         if ( out == v[i] )
  209.             error(E_INSITU,"zv_lincomb");
  210.  
  211.     out = zv_mlt(a[0],v[0],out);
  212.     for ( i = 1; i < n; i++ )
  213.     {
  214.         if ( ! v[i] )
  215.             error(E_NULL,"zv_lincomb");
  216.         if ( v[i]->dim != out->dim )
  217.             error(E_SIZES,"zv_lincomb");
  218.         out = zv_mltadd(out,v[i],a[i],out);
  219.     }
  220.  
  221.     return out;
  222. }
  223.  
  224.  
  225. #ifdef ANSI_C
  226.  
  227.  
  228. /* zv_linlist -- linear combinations taken from a list of arguments;
  229.    calling:
  230.       zv_linlist(out,v1,a1,v2,a2,...,vn,an,NULL);
  231.    where vi are vectors (ZVEC *) and ai are numbers (complex)
  232. */
  233.  
  234. ZVEC    *zv_linlist(ZVEC *out,ZVEC *v1,complex a1,...)
  235. {
  236.    va_list ap;
  237.    ZVEC *par;
  238.    complex a_par;
  239.  
  240.    if ( ! v1 )
  241.      return ZVNULL;
  242.    
  243.    va_start(ap, a1);
  244.    out = zv_mlt(a1,v1,out);
  245.    
  246.    while (par = va_arg(ap,ZVEC *)) {   /* NULL ends the list*/
  247.       a_par = va_arg(ap,complex);
  248.       if (a_par.re == 0.0 && a_par.im == 0.0) continue;
  249.       if ( out == par )        
  250.     error(E_INSITU,"zv_linlist");
  251.       if ( out->dim != par->dim )    
  252.     error(E_SIZES,"zv_linlist");
  253.  
  254.       if (a_par.re == 1.0 && a_par.im == 0.0)
  255.     out = zv_add(out,par,out);
  256.       else if (a_par.re == -1.0 && a_par.im == 0.0)
  257.     out = zv_sub(out,par,out);
  258.       else
  259.     out = zv_mltadd(out,par,a_par,out); 
  260.    } 
  261.    
  262.    va_end(ap);
  263.    return out;
  264. }
  265.  
  266.  
  267. #elif VARARGS
  268.  
  269. /* zv_linlist -- linear combinations taken from a list of arguments;
  270.    calling:
  271.       zv_linlist(out,v1,a1,v2,a2,...,vn,an,NULL);
  272.    where vi are vectors (ZVEC *) and ai are numbers (complex)
  273. */
  274. ZVEC  *zv_linlist(va_alist) va_dcl
  275. {
  276.    va_list ap;
  277.    ZVEC *par, *out;
  278.    complex a_par;
  279.  
  280.    va_start(ap);
  281.    out = va_arg(ap,ZVEC *);
  282.    par = va_arg(ap,ZVEC *);
  283.    if ( ! par ) {
  284.       va_end(ap);
  285.       return ZVNULL;
  286.    }
  287.    
  288.    a_par = va_arg(ap,complex);
  289.    out = zv_mlt(a_par,par,out);
  290.    
  291.    while (par = va_arg(ap,ZVEC *)) {   /* NULL ends the list*/
  292.       a_par = va_arg(ap,complex);
  293.       if (a_par.re == 0.0 && a_par.im == 0.0) continue;
  294.       if ( out == par )        
  295.     error(E_INSITU,"zv_linlist");
  296.       if ( out->dim != par->dim )    
  297.     error(E_SIZES,"zv_linlist");
  298.  
  299.       if (a_par.re == 1.0 && a_par.im == 0.0)
  300.     out = zv_add(out,par,out);
  301.       else if (a_par.re == -1.0 && a_par.im == 0.0)
  302.     out = zv_sub(out,par,out);
  303.       else
  304.     out = zv_mltadd(out,par,a_par,out); 
  305.    } 
  306.    
  307.    va_end(ap);
  308.    return out;
  309. }
  310.  
  311.  
  312. #endif
  313.  
  314.  
  315.  
  316. /* zv_star -- computes componentwise (Hadamard) product of x1 and x2
  317.     -- result out is returned */
  318. ZVEC    *zv_star(x1, x2, out)
  319. ZVEC    *x1, *x2, *out;
  320. {
  321.     int        i;
  322.     Real    t_re, t_im;
  323.  
  324.     if ( ! x1 || ! x2 )
  325.     error(E_NULL,"zv_star");
  326.     if ( x1->dim != x2->dim )
  327.     error(E_SIZES,"zv_star");
  328.     out = zv_resize(out,x1->dim);
  329.  
  330.     for ( i = 0; i < x1->dim; i++ )
  331.     {
  332.     /* out->ve[i] = x1->ve[i] * x2->ve[i]; */
  333.     t_re = x1->ve[i].re*x2->ve[i].re - x1->ve[i].im*x2->ve[i].im;
  334.     t_im = x1->ve[i].re*x2->ve[i].im + x1->ve[i].im*x2->ve[i].re;
  335.     out->ve[i].re = t_re;
  336.     out->ve[i].im = t_im;
  337.     }
  338.  
  339.     return out;
  340. }
  341.  
  342. /* zv_slash -- computes componentwise ratio of x2 and x1
  343.     -- out[i] = x2[i] / x1[i]
  344.     -- if x1[i] == 0 for some i, then raise E_SING error
  345.     -- result out is returned */
  346. ZVEC    *zv_slash(x1, x2, out)
  347. ZVEC    *x1, *x2, *out;
  348. {
  349.     int        i;
  350.     Real    r2, t_re, t_im;
  351.     complex    tmp;
  352.  
  353.     if ( ! x1 || ! x2 )
  354.     error(E_NULL,"zv_slash");
  355.     if ( x1->dim != x2->dim )
  356.     error(E_SIZES,"zv_slash");
  357.     out = zv_resize(out,x1->dim);
  358.  
  359.     for ( i = 0; i < x1->dim; i++ )
  360.     {
  361.     r2 = x1->ve[i].re*x1->ve[i].re + x1->ve[i].im*x1->ve[i].im;
  362.     if ( r2 == 0.0 )
  363.         error(E_SING,"zv_slash");
  364.     tmp.re =   x1->ve[i].re / r2;
  365.     tmp.im = - x1->ve[i].im / r2;
  366.     t_re = tmp.re*x2->ve[i].re - tmp.im*x2->ve[i].im;
  367.     t_im = tmp.re*x2->ve[i].im - tmp.im*x2->ve[i].re;
  368.     out->ve[i].re = t_re;
  369.     out->ve[i].im = t_im;
  370.     }
  371.  
  372.     return out;
  373. }
  374.  
  375. /* zv_sum -- returns sum of entries of a vector */
  376. complex    zv_sum(x)
  377. ZVEC    *x;
  378. {
  379.     int        i;
  380.     complex    sum;
  381.  
  382.     if ( ! x )
  383.     error(E_NULL,"zv_sum");
  384.  
  385.     sum.re = sum.im = 0.0;
  386.     for ( i = 0; i < x->dim; i++ )
  387.     {
  388.     sum.re += x->ve[i].re;
  389.     sum.im += x->ve[i].im;
  390.     }
  391.  
  392.     return sum;
  393. }
  394.  
  395. /* px_zvec -- permute vector */
  396. ZVEC    *px_zvec(px,vector,out)
  397. PERM    *px;
  398. ZVEC    *vector,*out;
  399. {
  400.     u_int    old_i, i, size, start;
  401.     complex    tmp;
  402.     
  403.     if ( px==PNULL || vector==ZVNULL )
  404.     error(E_NULL,"px_zvec");
  405.     if ( px->size > vector->dim )
  406.     error(E_SIZES,"px_zvec");
  407.     if ( out==ZVNULL || out->dim < vector->dim )
  408.     out = zv_resize(out,vector->dim);
  409.     
  410.     size = px->size;
  411.     if ( size == 0 )
  412.     return zv_copy(vector,out);
  413.     
  414.     if ( out != vector )
  415.     {
  416.     for ( i=0; i<size; i++ )
  417.         if ( px->pe[i] >= size )
  418.         error(E_BOUNDS,"px_vec");
  419.         else
  420.         out->ve[i] = vector->ve[px->pe[i]];
  421.     }
  422.     else
  423.     {    /* in situ algorithm */
  424.     start = 0;
  425.     while ( start < size )
  426.     {
  427.         old_i = start;
  428.         i = px->pe[old_i];
  429.         if ( i >= size )
  430.         {
  431.         start++;
  432.         continue;
  433.         }
  434.         tmp = vector->ve[start];
  435.         while ( TRUE )
  436.         {
  437.         vector->ve[old_i] = vector->ve[i];
  438.         px->pe[old_i] = i+size;
  439.         old_i = i;
  440.         i = px->pe[old_i];
  441.         if ( i >= size )
  442.             break;
  443.         if ( i == start )
  444.         {
  445.             vector->ve[old_i] = tmp;
  446.             px->pe[old_i] = i+size;
  447.             break;
  448.         }
  449.         }
  450.         start++;
  451.     }
  452.     
  453.     for ( i = 0; i < size; i++ )
  454.         if ( px->pe[i] < size )
  455.         error(E_BOUNDS,"px_vec");
  456.         else
  457.         px->pe[i] = px->pe[i]-size;
  458.     }
  459.     
  460.     return out;
  461. }
  462.  
  463. /* pxinv_zvec -- apply the inverse of px to x, returning the result in out
  464.         -- may NOT be in situ */
  465. ZVEC    *pxinv_zvec(px,x,out)
  466. PERM    *px;
  467. ZVEC    *x, *out;
  468. {
  469.     u_int    i, size;
  470.     
  471.     if ( ! px || ! x )
  472.     error(E_NULL,"pxinv_zvec");
  473.     if ( px->size > x->dim )
  474.     error(E_SIZES,"pxinv_zvec");
  475.     if ( ! out || out->dim < x->dim )
  476.     out = zv_resize(out,x->dim);
  477.     
  478.     size = px->size;
  479.     if ( size == 0 )
  480.     return zv_copy(x,out);
  481.     if ( out != x )
  482.     {
  483.     for ( i=0; i<size; i++ )
  484.         if ( px->pe[i] >= size )
  485.         error(E_BOUNDS,"pxinv_vec");
  486.         else
  487.         out->ve[px->pe[i]] = x->ve[i];
  488.     }
  489.     else
  490.     {    /* in situ algorithm --- cheat's way out */
  491.     px_inv(px,px);
  492.     px_zvec(px,x,out);
  493.     px_inv(px,px);
  494.     }
  495.     
  496.     
  497.     return out;
  498. }
  499.  
  500. /* zv_rand -- randomise a complex vector; uniform in [0,1)+[0,1)*i */
  501. ZVEC    *zv_rand(x)
  502. ZVEC    *x;
  503. {
  504.     if ( ! x )
  505.     error(E_NULL,"zv_rand");
  506.  
  507.     mrandlist((Real *)(x->ve),2*x->dim);
  508.  
  509.     return x;
  510. }
  511.